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Abstract 

We have performed a high-precision Monte Carlo study of the dynamic critical be- 
havior of the Swendsen-Wang algorithm for the three-dimensional Ising model at the 
critical point. For the dynamic critical exponents associated to the integrated auto- 
correlation times of the "energy-like" observables, we find z-^j^f = Zi n t,£ = £int,£' 
0.1")!) ± 0.005 ± 0.025, where the first error bar represents statistical error (68% con- 
fidence interval) and the second error bar represents possible systematic error due to 
corrections to scaling (68% subjective confidence interval). For the "susceptibility-like" 
observables, We find ^intjAI 2 — ^'mt,S2 — 

0.443 ± 0.005 ± 0.030. For the dynamic criti- 
cal exponent associated to the exponential autocorrelation time, we find. .s CX p ps 0.481. 
Our data are consistent with the Coddington-Baillie conjecture zsw = P/v ~ 0.5183, 
especially if it is interpreted as referring to z exp . 



PACS codes: 05.50.+q, 05.10.Ln, 64.60. Cn, 64.60.Ht. 



Key Words: Ising model; Potts model; Swendsen-Wang algorithm; cluster algorithm; 
Monte Carlo; autocorrelation time; dynamic critical exponent. 



1 Introduction 



Monte Carlo (MC) simulations [1-7] have become a standard and powerful tool for gain- 
ing nonperturbative insights into statistical-mechanical systems and lattice field theories. 
However, their practical success is severely limited by critical slowing-down: the autocor- 
relation time t — that is, roughly speaking, the time needed to produce one "statistically 
independent" configuration — diverges near a critical point. More precisely, for a finite 
system of linear size L at criticality, we expect a behavior r ~ L z for large L. The power 
z is a dynamic critical exponent, and it depends on both the system and the Monte Carlo 
algorithm. 

Local Monte Carlo algorithms (such as single-site Metropolis or heat bath) generally have 
a dynamic critical exponent z > 2. This makes it very hard to get high-precision data very 
close to the critical point on large lattices. 

In some cases, a much better dynamical behavior can be obtained by including non-local 
moves, such as cluster flips. 1 In particular, the Swendsen-Wang (SW) cluster algorithm [9] 
for the ferromagnetic g-state Potts model achieves a significant reduction in z compared 
to the local algorithms: one has z between and 1, where the exact value depends on q 
and on the dimensionality of the lattice. The most favorable case is the two-dimensional 
(2D) Ising model (q = 2), for which the best currently available numerical estimate is 
z = 0.222 ±0.007 [10] (see also the discussion below). In other cases, the performance of the 
SW algorithm is somewhat less impressive but still quite good: e.g., z = 0.514 ± 0.006 for 
the 2D 3-state Potts model [11], z « 1 for the 2D 4-state Potts model [12,13], z « 0.45-0.75 
for the 3D Ising model [9,14-20], and z » 0.86-1 for the 4D Ising model [16,21-24]. Clearly, 
we would like to understand why the SW algorithm works so well in some cases and less well 
in others. We hope in this way to obtain new insights into the dynamics of non-local Monte 
Carlo algorithms, with the ultimate aim of devising new and more efficient algorithms. 

There is at present no adequate theory for predicting the dynamic critical behavior of an 
SW-type algorithm. However, Li and Sokal [25] have proven that the autocorrelation times 
of the Swendsen-Wang algorithm for the g-state Potts ferromagnet are bounded below by a 
multiple of the specific heat: 

Tmt,M, 7int,£, Ti nt ,£', T exp > Const X C H ■ (1.1) 

Here M is the bond occupation in the SW algorithm, S is the energy, £' is the nearest- 
neighbor connectivity, and Ch is the specific heat; r int and r exp denote the integrated and 
exponential autocorrelation times, respectively [4]. As a result one has for the dynamic 
critical exponents 

z mt,Ni z mt,£, z xat,£'i z exp > Oi / V , (1.2) 

where a and v are the standard static critical exponents. Thus, the SW algorithm cannot 
completely eliminate the critical slowing-down if the specific heat is divergent at criticality. 2 

1 See [4,8] for reviews of collective- mode Monte Carlo methods. 

2 The bound / (|1.2|) has also been proven to hold [12] for the direct SW-type algorithm [26] for the 
Ashkin- Teller (AT) model [27,28], which interpolates between the 2-state (Ising) and 4-state Potts models. 
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The physical mechanism underlying the Li-Sokal proof is the slow evolution of the bond 
occupation Af: its mean-square change per SW iteration is of order V (volume), while its 
variance in the equilibrium distribution is of order VCh] this leads to (jl.ljl . In addition, 
Salas and Sokal [11] have proven, under a very mild and eminently plausible condition 3 , that 
all three "energy-like" observables have the same dynamic critical exponent: 

z int,M — z int,£ = z mt,S' ■ (1-3) 

For details on all these results, see [11, Section 2.2]. 

One important question is whether the Li-Sokal bound (|1.1|) / (|1.2j) is sharp or not. An 
affirmative answer would imply that we could use the bound to predict the dynamic crit- 
ical exponent (s) z given only the static critical exponents of the system. There are three 
possibilities: 

i) The bound (|1.1|) is sharp (i.e., the ratio t/Ch is bounded), so that (jl.2|) is an equality. 

ii) The bound is sharp modulo a logarithm (i.e., t/Ch ~ log p L for p > 0). 

iii) The bound is not sharp (i.e., t/Ch ~ L p for p > 0), so that (jl.2|) is a strict inequality. 

The best currently available data concern the two-dimensional Potts models with q — 2, 3, 4, 
for which we have 

q = 2: a/z/ = 0(xlog) z- mt £' = 0.222 ± 0.007 [10] 
q = 3: a/v = 2/5 z- m t,e> = 0.514 ± 0.006 [11] 

q = 4: a/v = 1 (x log~ 3/2 ) z int , £ = 0.876 ± 0.011 [13] 

Here the values of a/v are exact [13,29,30], while the values of z are the best available 
numerical estimates from pure power-law fits. Note, however, that the estimate of z for 
q = 4 cannot be correct, as it violates the Li-Sokal bound (|1.2j) : presumably it is corrupted 
by the same multiplicative logarithmic corrections that afflict the specific heat. 4 For this 
reason, the papers [10,11,13] analyzed also the ratio t/Ch in order to test directly the 
sharpness of the Li-Sokal bound. It was found that the data for q = 2,3,4 are consistent 
with two scenarios: either the Li-Sokal bound is non-sharp by a very small power (p w 0.06- 
0.12), or else it is sharp modulo a logarithm (possibly with power p — 1). Not surprisingly, 
it is exceedingly difficult to distinguish numerically between these two scenarios. 

For the three- and four- dimensional Ising models, by contrast, the Li-Sokal bound (jl.2j) 
is clearly not sharp: the numerical estimates of z are much larger than a/v. It follows 
that another physical mechanism, beyond the one captured in the Li-Sokal proof, must be 
principally responsible for the critical slowing-down in these cases; but it is far from clear 
what this mechanism is. A natural first step towards identifying this mechanism would be to 
obtain accurate numerical estimates for z in the three- and four-dimensional Ising models. 

3 The condition is that the normalized autocorrelation functions of the bond occupation at time lags 1 
and 2, i.e. ptva/^I) and pj\f^f(2), should be bounded away from zero as the lattice size L tends to infinity. 

4 Indeed, a pure power-law fit to the Monte Carlo data for the specific heat yielded ajv = 0.770±0.008 [13] . 
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Unfortunately, the numerical results in the four- dimensional case are almost nonexistent 5 , 
and in the three-dimensional case are wildly contradictory. In chronological order, the results 
for the three-dimensional Ising model are as follows: 

1987: -2 exp ,£- = 0.75 ± 0.01, based on unspecified lattice sizes [9] 
1989: z- mt ,e = z intM 2 = 0.50 ± 0.03, based on L < 64 [14] 

1990: ;z cxp a = 0.46 ± 0.01 (Ci = size of the largest cluster), based on L < 89 
extrapolated to L = oo [15] 6 

1992: z intt£ = 0.54 ± 0.02, based on L < 48 [16] 

1993: z exp = 0.61 ± 0.02, z int ,s = 0.58 ± 0.02, z intM 2 = 0.57 ± 0.02 and z int> \ M \ = 
0.55 ± 0.02, based on L < 32 [18] 

1993: z iati£ = 0.49 ± 0.02 and z intt]Ml = 0.48 ± 0.01, based on L < 60 [19] 

2002: z- mt>£ « 0.50, based on L < 128 [20] 7 

The reason for the discrepancies is unclear, but there does seem to be a tendency for the esti- 
mates of z to decrease as larger lattices are used — an effect that could easily be understood 
as arising from corrections to scaling. 

The purpose of this paper is to restudy the dynamic critical behavior of the SW algorithm 
for the three-dimensional Ising model, using much larger lattices (up to L — 256), vastly 
higher statistics (well over 10 7 SW iterations at each lattice size), and a careful finite-size- 
scaling analysis. For the dynamic critical exponents associated to the integrated autocorre- 
lation times of the "energy-like" observables, we find 

Zint,N = Zint,£ = = 0.459 ± 0.005 ± 0.025 , (1.4) 

where the first error bar represents statistical error (68% confidence interval) and the second 
error bar represents possible systematic error due to corrections to scaling (68% subjective 

5 The only numerical study of which we are aware is [16], which yielded z- mt .£ = 0.86 ± 0.02, based on 
lattices of size 4 < L < 16 (which by present-day standards are much too small). In addition, there have been 
numerical [22, 23] and analytic [22-24] studies of the Swendsen-Wang algorithm for the Ising ferromagnet 
on the complete graph (also known as the Curie- Weiss or "mean- field" model), which indicate z = 1. This 
model is 'presumed to lie in the same dynamic universality class as the Ising model on a regular lattice of 
dimension d > 4, but high-precision numerical tests of this quite plausible conjecture are lacking. 

6 Wang [15] also studied the magnetization (including sign) of the largest cluster in a variant of the 
Swendsen-Wang algorithm in which the largest cluster is not flipped. But it is far from clear whether 
this observable corresponds to any observable in the standard Swendsen-Wang algorithm. Indeed, Wang 
found that the exponential autocorrelation time of this observable is several times larger than that of C±. 
And since there is good reason to believe (see Section 15.21 below) that C\ does indeed have a significant 
overlap with the slowest mode in the Swendsen-Wang algorithm, this suggests that Wang's observable is not 
interpretable within the standard Swendsen-Wang algorithm, but rather represents a new slow mode in the 
variant algorithm. 

pure power-law fit to the raw data of Wang, Kozan and Swendsen [20] yields a decent \ 2 if (and only 
if) imin > 32. Our preferred fit is L min = 32, and yields z int , £ = 0.502 ± 0.012 ( X 2 = 0.440, 1 DF, level = 
50.7%). We thank Jian-Sheng Wang for supplying us with these raw data. 
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confidence interval). For the "susceptibility-like" observables, we find 

z-mt,M 2 = z mt,s 2 = 0.443 ± 0.005 ± 0.030 . (1.5) 

Finally, for the dynamic critical exponent associated to the exponential autocorrelation time, 
we obtain the rough estimate 

z exp « 0.481 . (1.6) 

It is possible that some or all of these exponents are in fact exactly equal. 

The present paper is organized as follows: Section previews the basics of the Swendsen- 
Wang algorithm and the definitions of autocorrelation times and observables. In Sectional 
we discuss our methods of statistical data analysis. In Section 0] we summarize our Monte 
Carlo simulations. In Section |5] we present the analysis of our dynamic data. In Section El 
we discuss our results in the light of various conjectures that have been made by previous 
workers. The static data from our simulations will be analyzed in a separate paper [31]. 



2 Basic set-up and notation 



2.1 Potts model and Swendsen— Wang algorithm 

The g-state Potts model assigns to each lattice site i a spin variable cr, taking values in 
the set {1,2, ... ,q}; these spins interact through the reduced Hamiltonian 



^Potts = — p y^(^i,(Tj ~ i) 

(ij) 



(2.1) 



where the sum runs over all the nearest-neighbor pairs (ij) (each pair counted once). To 
simplify the notation we shall henceforth write 5 aijCr . = 5 ab for a bond b = (ij). The 
ferromagnetic case corresponds to f3 > 0. The partition function is defined as 



-H 



£exp 

M 



i) 



(2.2) 



Finally, the Boltzmann weight of a configuration {a} is given by 



WWW) 



exp 



0£(^ 



(2.3) 



where p — 1 — e _/3 . 

The idea behind the Swendsen- Wang (SW) algorithm [9, 32] is to decompose the Boltz- 
mann weight ()2.3|) by introducing new dynamical variables n& = 0, 1 living on the bonds of 
the lattice, and to simulate the joint model of old and new variables by alternately updating 
one set of variables conditional on the other set. The Boltzmann weight of the joint model 
is 

W ioint ({a}, {n}) = I J] [(1 - p)S nb>0 + P 5 ai S nb>1 ) . (2.4) 
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The marginal distribution of (|2.4|) with respect to the spin variables reproduces the Potts- 
model Boltzmann weight ()2.3j) . The marginal distribution of ()2.4|) with respect to the bond 
variables is the Fortuin-Kasteleyn [33-35] random-cluster model with parameter q: 



Wrc({™}) = ^ 



n 



P 



n 

b: nj,=0 



C({n}) 



(2.5) 



where C({n}) is the number of connected components (including one-site components) in 
the graph whose edges are the bonds with n& = 1. 

We can also consider the conditional probabilities of the joint distribution ([2.4)1 . The 
conditional distribution of the {n} given the {a} is as follows: independently for each bond 
b = (ij), one sets n b = when a.i ^ Oj, and sets n b = and 1 with probabilities 1 — p and 
p when Oi = a,y Finally, the conditional distribution of the {a} given the {n} is as follows: 
independently for each connected cluster, one sets all the spins Oi in that cluster equal to 
the same value, chosen with uniform probability from the set {1,2, ... ,q}. 

The Swendsen-Wang algorithm simulates the joint probability distribution ()2.4|) by al- 
ternately applying the two conditional distributions just described. That is, we first erase 
the current {n} configuration, and generate a new {n} configuration from the conditional 
distribution given {a}; we then erase the current {a} configuration, and generate a new {a} 
configuration from the conditional distribution given {n}. A single step of the SW algorithm 
consists of these two "half-steps" . 



2.2 Autocorrelation functions and autocorrelation times 

Let O be any observable (i.e. any function of {a} and {n}), and let 0(t) be its evolution 
in Monte Carlo time (where one unit of time corresponds to a single step of the Swendsen- 
Wang algorithm). The unnormalized autocorrelation function associated to the observable 
O is defined as 

Coo(t) = (0(s)0(s + t))-(0) 2 , (2.6) 

where the expectations are taken in equilibrium. The corresponding normalized autocorre- 
lation function is defined as 

( x _ C o{t) C 00 (t) 

Poo{t = 7^ — 77^- - • (2- 'J 

Coo(0) var(C') 

The integrated and exponential autocorrelation times associated to the observable O are 
defined as 8 



Tint,© 



7> E Pooit) (2.8) 



2 

t 



8 For a general Markov chain, the "lim" in l|2.9() should strictly speaking be replaced by "limsup", and 
poo(t) should be replaced by its absolute value. But in the Swendsen-Wang algorithm it can be proven that 
the limit really exists, and that poo(t) > for all t; this follows from the spectral representation [25] [11, 
Section 2.2]. 
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r exp ,o = ,lim — — — — (2.9) 



\t 

\t\-^o log poo (t) 

Finally, the exponential autocorrelation time of the system is defined as 



r exp = sup r exPi o , (2.10) 
o 

where the supremum is taken over all observables O. This autocorrelation time measures 
the decay rate of the slowest mode of the system. All observables that are not orthogonal to 
this slowest mode satisfy r exP) e> = r cxp . 

It is important to remember that there is not just one autocorrelation time, but many: 
namely, r cxp as well as r- mt) o for each observable O. In all but the most trivial Markov chains, 
these autocorrelation times are not equal. Correspondingly, there are many dynamic critical 
exponents: namely, z exp as well as Zi n t,e> for each observable O. These exponents may in some 
cases be equal (i.e., the corresponding autocorrelation times may scale proportionally as the 
critical point is approached), but they need not be; this is a detailed dynamical question, 
and the answer will vary from model to model. 

2.3 Observables to be measured 

As just explained, the Swendsen-Wang algorithm is most naturally defined in the general 
context of the g-state Potts ferromagnet. It is therefore most convenient and natural to use 
a formalism that is valid for arbitrary q; at the end we can specialize to the Ising case q = 2. 

The nicest "geometric" representation of Potts spins is the hypertetrahedral represen- 
tation, defined as follows: Let {e^} 9 a=1 be unit vectors in IR 9-1 satisfying • = 
(qS a/3 — l)/(<7 — 1). Geometrically, these vectors point from the center to the vertices of a 
unit hypertetrahedron in IR 9-1 . We then represent a Potts spin o~ x G {1,2, ... ,q} by the 
unit vector cr x = in M 9_1 . This representation captures the S q (permutation group) 
symmetry of the Potts Hamiltonian (|2.1|) . and for q = 2 it reduces to the usual representation 
of Ising spins cr x — ±1. We have in particular 

<r*-"v= q5 ^~ X , (2.11) 
so that the Potts Hamiltonian can be written equivalently as 

H = -PpoUb y, frr f ,(r.j + const (2.12a) 

= -Aetr y, <Ti • cfj + const (2.12b) 

(ij) 

where 

Aetr — A'otts • (2-13) 

Q 

For q = 2 this yields Asing = A'otts /2, where Asing = Aetr corresponds to the usual Ising 
normalization for the inverse temperature. 
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Let us now consider the g-state Potts ferromagnet on a (/-dimensional periodic hypercubic 
lattice of linear size L. We write V = L d for the number of sites, and B = dL d for the number 
of bonds. We shall consider the following observables: 

• (minus) the total energy 

S = ^(T x -(Ty (2.14a) 

= j2 q6ax, 2^ 1 ( 2 - 14b ) 

{xy) ^ 

where the sum runs over all the nearest-neighbor pairs (xy) (each pair counted once). 

• the bond occupation 

A/" = J2 n *y ( 2 - 15 ) 

• the nearest-neighbor connectivity (which is an energy-like observable [11]) 

{xy) 

where 7^ equals 1 if both ends of the bond (xy) belong to the same cluster, and 
otherwise. More generally, the connectivity 7^ can be defined for an arbitrary pair i,j 
of sites: 

/r\\ f 1 if z is connected to j , , 

li3\\ n i) I jf i i s no t connected to j 

We shall also use higher connectivities, such as 

, r -1 \ f 1 if i, j.k, I are all connected together , n , oX 

^"« n » = 1 otherwise ( 2 ' 18 ) 

• the squared magnetization 

M 2 = (j>^ ( 2 - 19a ) 



X 



ct=l \ x 



• higher powers of the magnetization 

M 2n = {M 2 ) n (2.20) 
(In this paper, we measured only M 2 and M A .) 
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the square of the Fourier transform of the spin variable at the smallest allowed non-zero 
momentum 



1 d 



d 



Q 



d q 



3 = 1 ce=l 



E^ 



(2.21a) 
(2.21b) 



where (x±, x 2 , ■ ■ ■ , xj) are the Cartesian coordinates of point x. Note that T is normal- 
ized to be comparable to its zero- momentum analogue M. 2 . 

• the number of clusters (= connected components) and the mean-square and mean- 
fourth-power size of the clusters 



C = S 



c 

C x,y 

E#( C ') 4 = E 7*y«» 



(2.22) 
(2.23) 
(2.24) 



x,y,u,v 



where the sum runs over all the clusters C of activated bonds, and #(C) is the number 
of sites in the cluster C. 

• the size C« of the ith largest cluster (C\ > C 2 > C 3 > . . .). In this work we measured 
only Ci, C 2 and C 3 . 

From these observables we compute the following expectation values: 

• (minus) the energy density E per bond 

1 



e = 5 <£> , 



(2.25) 



where B = dV is the number of bonds in the lattice, so that a perfectly disordered 
(resp. ferromagnetically ordered) state has E = (resp. E = 1) 



• the specific heat per bond 

• the magnetic susceptibility 



^var(£) = l[<£ 2 >-<£> 2 ] 



(2.26) 



(2.27) 
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the correlation function at momentum (27r/L, 0, . . . , 0) 



F = i(J-) (2.28) 



* = o^Lm (I- 1 )' 72 (2-29) 



the second-moment correlation length 

2sm(7r/L) \F 

• the mean number of clusters 

So = (S ) (2.30) 

• the mean size of the ith largest cluster 

a = (Ci) (2.31) 

For each observable O discussed above, we have measured its autocorrelation function 
Pooif) and have used this to estimate the corresponding integrated autocorrelation time 
Ti ntj e). In Section El we explain how we derived estimates of the mean values and the error 
bars for both static and dynamic quantities. 

Remarks. 1. Using the Fortuin-Kasteleyn identities [4,33-35], which arise from the for- 
mulae for conditional expectations in the joint measure ()2.4|) . it is not difficult to show 
that 

(n xy ) = p(5 ax , ay ) (2.32) 
(<T x -cr y ) = ( lxy ) (2.33) 



-Zlxyuv)) (2.34) 



and hence that 



Q 

(£) = (£') (2.36) 
(M 2 ) = (S 2 ) (2.37) 

(M 4 ) = q -^(Sl) - -^-<S 4 > (2.38) 
q — 1 q — 1 

where p = 1 — e _/3potts is the Swendsen-Wang bond probability and B = dL d is the number 
of bonds in the lattice. As a check on the correctness of our simulations, we have tested 
these identities to high precision, in the following way: Instead of comparing directly the left 
and right sides of each equation, which are strongly positively correlated in the Monte Carlo 
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simulation, a more sensitive test is to define new observables corresponding to the differences 
(i.e., £ — £' and so forth). Each such observable should have mean zero, and the error bars 
on the sample mean can be estimated using the standard error analysis outlined in Sectional 
below. The comparison to zero yields the following x 2 values: 



For (j2~35J) 
For (mj) 
For (I2~3T|) 
For (EH| 



X 2 = 160.95 (182 DF, level = 87%) 
X 2 = 152.63 (182 DF, level = 94%) 
X 2 = 188.28 (182 DF, level = 36%) 
X 2 = 190.47 (182 DF, level = 32%) 



(2.39) 
(2.40) 
(2.41) 
(2.42) 



Here we have treated each independent run (see Section HJ) as a separate data point; DF 
means the number of degrees of freedom (i.e. the number of independent runs), and "level" 
means the confidence level of the fit (defined at the beginning of Section below). The 
agreement is excellent. 

2. As a further check on the correctness of our simulations, we have computed both sides 
of the identity 



Paw(!) 



1 - 



q 2 (l-p)E 



(2.43) 



(q - l) 2 V C H + q 2 (l-p)E 

proven in [25, equation 7] and [11, equation (2. 16)]). 9 This is a highly nontrivial test, as it 
relates static quantities (energy and specific heat) to a dynamic quantity (autocorrelation 
function of the bond occupation at time lag 1). We have also checked with great accuracy 
the identities [11] 



C ee {t) 
Pes{t) 

P£'£'{t) 



g-1 

Pa/-a/-(1) 

C se (t+1) 

Peeit + l) 
Pes(l) 



(2.44) 

(2.45) 
(2.46) 
(2.47) 



3 Statistical methods 

In this paper we are aiming at extremely high precision for both static and dynamic 
quantities; and furthermore we need to disentangle the effects of statistical errors from the 
effects of systematic errors due to corrections to scaling. For this, it is essential to obtain 
accurate estimates not only of the static and dynamic quantities of interest, but also of their 
error bars: in this way we will be able (see Section |SJ) to perform x 2 tests that provide an 
objective measure of the goodness of fit in each scaling Ansatz. 

9 Unfortunately, three different normalizations of the energy are used in [25], in [11], and in the present 
paper. 
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In this section we review briefly how we performed the statistical analysis of our raw 
Monte Carlo data. In particular, we describe how to compute the estimators for the mean 
value and the variance of both static and dynamic quantities. These methods are based on 
well-known results of time-series analysis [36,37]. More details on the methods used here 
can be found in [38, Appendix C], [4, Section 3] and [11, Section 4]. 

Let us consider a generic observable O, whose mean is equal to po- Its corresponding un- 
normalized and normalized autocorrelation functions are denoted by Coof) = (^(0)^00) — 
(O) 2 and poo(t) = Coo(t)/Coo(fy, respectively. We also define the integrated autocorrela- 
tion time 

j oo 

Ti„ t ,o = Poo{t) . (3.1) 

t=— oo 

Given a sequence of n Monte Carlo measurements of the observable O — call them 
. . . , O n } — the natural estimator of the mean \i is the sample mean 



1 n 

O = -VOi. (3.2) 



n . 



This estimator is unbiased and has a variance 

1 n 

var(O) = -V(7 00 (r- S ) (3.3a) 



r,s=l 

1 ^ I > 



E (l--)Coo(t) (3.3b) 



n * — ' \ n 

t=-(n-l) 

~ - 2r intjC , Co O (0) for n > r intjC , (3.3c) 

This means that the variance is a factor 2ri nti o larger than it would be if the measurements 
were uncorrelated. It is, therefore, very important to estimate the autocorrelation time r inti e> 
in order to ensure a correct determination of the error bar on the (static) quantity fio- 
The natural estimator for the unnormalized autocorrelation function Cooif) is 

n-\t\ 



if the mean ji is known, and 



1 1 i=l 



if the mean \iq is unknown. We emphasize that, for each t, the estimators C 00(f) and 
Coo{t) are random variables [in contrast to Coo(t), which is a number]. The estimator 
C 00(f) is unbiased, and C 00(f) is biased by terms of order 1/n. The covariance matrices of 
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Cqo an d Coo are the same to leading order in the large-n limit (i.e., n ^> r int C )), and we 
have [36,37] 

cov {Coo (t), Coo («)) = - [Coo (m) C o {m + u - t) + C 00 {m + u) C o {m - t) 



n 

i=—oo 



+K,(t, m, m + u)] + o ( — J , (3.6) 

W 

where t, u > and k is the connected 4-point autocorrelation function 

n(r,s,t) = ((Oi- fio)(O i+r - fi )(O i+s - fi )(O i+t - fi )) 
-C o(r)C o(t ~ s) - Coo(s)C o(t ~ r) 
-C o(t)Coo(s - r) . (3.7) 

The natural estimator for the normalized autocorrelation function poo{t) is 



if the mean p is known, and 



C oo (0) 



if the mean p a is unknown. The estimators poo(t) and p oi^) are biased by terms of order 
1/n, as a result of the ratios of random variables in ()3.8|) / ()3.9|) . The covariance matrices of 
'poo an d f>oo are the same to leading order in 1/n. If the process is Gaussian, this covariance 
matrix is given in the large-n limit by [37] 

- oo 

cov(p 00 (t), poo{u)) = - } j [poo(m)poo(m + t-u) + poo(m + u)p o(m - t) 

m=— oo 

+2poo(t)poo(u)p 2 00 (m) - 2poo(t)poo(m)poo(m - u) 
-2poo(u)poo(m)poo(m -t)]+o ( - J (3.10) 



?? 



for i, u > 0. If the process is not Gaussian, then there are additional terms proportional to 
the fourth cumulant K,(m,t,t — u). The simplest assumption is to consider the stochastic 
process to be "not too far from Gaussian", and drop all the terms involving k. If this 
assumption is not justified, then we are introducing a bias in the estimate of this covariance. 
Finally, we shall take the estimator for the integrated autocorrelation time to be [38] 



1 M 

Tint,© = g Yl POo{t) (3.11) 



t=-M 
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[or the same thing with Poo(t)] where M is a suitably chosen number. The reason behind 
the cutoff M is the following: if we were to make the "obvious" choice M = n + 1, then the 
resulting estimator would have a variance of order 1 even in the limit n — > oo; this is because 
the terms f>oo{t) with large t have errors (of order 1/n) that do not vanish as t grows [cf. 
f)3.1Ujl ]. and their number is also large (~ n). Taking M <C n restores the good behavior of 
the estimator as n — > oo. The bias introduced by this rectangular cutoff 10 is given by 



bias(f inti o) = -- ^2 Pooit) + o (- J 

\t\>M ^ ' 



(3.12) 



The variance of the estimator T mX ,,o can be computed from the covariance ()3.10|) ; the final 
result is [38] 

2(2M+1) 2 

var(r inti0 ) « r int , (3.13) 

where the approximation r intj e) <C M -C n has been made. A good (self-consistent) choice of 
M is the following [38]: let M be the smallest integer such that M > CT^pt o(M), where c is a 
suitable constant. If the normalized autocorrelation function is roughly a pure exponential 11 , 
then a choice in the range c ~ 6-8 is reasonable. Indeed, if we take poo{t) — e~ l ^ T and 
minimize the mean-square error 

MSE(? int ,c>) = bias(f intiC ,) 2 + var(f int)0 ) (3.14) 

using ()3.12|) / ()3.13|) , we find that the optimal window width is 

For n/r ~ 10 8 (resp. 10 6 , 10 4 ), we have M opt /r w 8.86 (resp. 6.56, 4.26). In this paper we 
used c = 8 for the observables J\f, £, £', Ai 2 , S 2 and C±, whose autocorrelation functions are 
close to a pure exponential (see Section IB^j) : c = 10 for So] and c = 15 for C 2 and C3. 

As noted above, we expect the estimator T inti e> to have a bias of order T inti e>/n, due to the 
nonlinearities in (|3.8J) / (J3.9|) . 12 To make this bias negligible we need long runs. It has been 
shown empirically that this procedure works fairly well when n > 10 4 r inti £) [4]. 

Remarks. 1. For the specific heat Cjj arid the correlation length £, which are "composite" 
quantities (i.e. not merely the mean value of a single observable), the estimation of the error 
bars is a bit more complicated. One method is described in [11, Section 4]; a slightly better 
method, based on the analysis of the cross-correlation matrix, is described in [31]. 

10 We could use more general cutoff functions, but this rectangular cutoff is the most convenient for the 
present purposes. 

11 This is confirmed here for the observables AT, £, Ai 2 , S2 and C\ (see Section 15121 . Similar behavior is 
found in the SW algorithm for the two-dimensional Ising [10], 3-state Potts [11] and 4-state Potts models [12]. 

12 The bias on the estimator Ti nti e> also induces a bias on the estimated variance (|3.3(l of the sample mean 
O. This bias is of order 1/n 2 , i.e. a factor 1/n down from the variance (|3.3|l itself. 
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2. On most lattices we made a number of independent runs, rather than one long run 
(see Section |U below). Our best estimate of each autocorrelation function pooif) was then 
obtained by averaging the estimates f>oo(t) fr° m the individual runs, with weights propor- 
tional to the run lengths. Finally, the windowing procedure was performed on the resulting 
best estimate of poo{t)- This is a better procedure than performing the windowing on each 
run separately. 

3. As a check on the correctness of the error bars produced by our time-series-analysis 
method, we also used an analysis method based on independent "bunches" [11, Section 4.2]. 
The error bars ranged from « 50% to ~ 115% of those produced by the time-series-analysis 
method, averaging around 80%. The fluctuations are not surprising, as the bunch-method 
error bar has a statistical fluctuation of order l/y/m, where m is the number of bunches (in 
our case ranging from 10 to 35). However, the systematic tendency toward smaller error bars 
suggests that our time-series-analysis method may be slightly overestimating the error bars, 
probably due to neglect of the non-Gaussian terms in ()3.10|) . Therefore, the true statistical 
error bars on our raw data and on our exponent estimates may be slightly smaller than those 
reported in this paper. This issue deserves a more detailed investigation in the future. 



4 Description of the simulations 

We implemented the Swendsen-Wang algorithm for the nearest-neighbor three-dimensional 
Ising model on an L x L x L simple-cubic lattice with periodic boundary conditions. We 
performed all our runs at /?i sing = 0.22165459 (i.e. /3p otts = 0.44330918), which is Blote et 
a/.'s [39] best estimate of the critical temperature and is very near to the estimates by other 
workers [40,41] (see also the review [42]). We studied lattice sizes L — 4, 6, 8, 12, 16, 24, 
32, 48, 64, 96, 128, 192, 256 and performed between 2.9 x 10 7 and 5 x 10 8 SW iterations 
for each lattice size (see Table C]). The total data set at each L corresponds to ~ 10 6 r on 
the largest lattices (L = 192,256), at least 10 7 r on all L < 64, and nearly 10 8 r at L — 16 
(see again Table In all cases, the statistics are high enough to permit a high accuracy in 
our estimates of the static (error ~ 0.01-0.12%) and dynamic (error ~ 0.1-0.5%) quantities. 
Our results for the principal static observables are reported in Table El and for the dynamic 
quantities in Tables El and |U 

The initial configuration of each run was either random or ordered, and we discarded 
the first 10 5 iterations from each run in order to allow the system to reach equilibrium; this 
discard interval is in all cases greater than 4000 T inti £', which is more than sufficient. 13 We 

13 Such a discard interval might seem to be much larger than necessary: 100r eX p would usually be more 
than enough. However, there is always the danger that the longest autocorrelation time in the system (r eX p) 
may be much larger than the longest autocorrelation time that one has measured, because one has failed 
to measure an observable having sufficiently strong overlap with the slowest mode. As an undoubtedly 
overly conservative precaution against the possible (but unlikely) existence of such a (vastly) slower mode, 
we decided to discard 10 5 iterations. In most cases this amounts to less than 10% of the run, thus reducing 
the accuracy on our final estimates by less than 5%. Unless there exists a vastly slower mode of which we 
are unaware, our data yield T- mt ^' /r cxp ss 0.9-1 for this algorithm (see Section and Table [S] below) . So 
the discard interval is greater than 4000r oxp . 
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checked that random and ordered initial conditions gave identical results, within statistical 
error. On some of the smaller lattices, we made a single long run of 10 8 iterations; on other 
lattices, we averaged the data from several (anywhere from 2 to 46) individual runs of at least 
10 6 iterations each (except for a small number of runs of length 5 x 10 5 at L = 96, 192, 256). 
In all cases we discarded the first 10 5 iterations of each run. The individual runs (minus 
the discard) are all of length greater than 20000ri nti £', which is long enough to allow a good 
determination of the dynamic quantities. 

Our program was written in Fortran 77 and run on a 1266 MHz Pentium III Tualatin 
processor using the g77 Fortran compiler. Our program requires approximately 42L 3 bytes 
memory. The CPU time required by our program ranges from 0.39 to 0.90 L 3 /xs/iteration, 
depending on the lattice size (see Table Q). The sharp rise in CPU time per spin on very 
small lattices arises from the "fixed costs" of the algorithm (i.e. those that do not scale with 
the volume). The slow rise in CPU time per spin on larger lattices arises from the "cache 
misses" that occur, due to the nonlocal nature of the Swendsen-Wang algorithm, when the 
lattice no longer fits in the 512 KB cache. The total CPU time used in these runs was 
approximately 17.8 years. 

In the first version of our program, the random numbers were supplied by a linear con- 
gruential generator 

x n+ i = ax n + c (mod m) (4.1) 

with modulus m = 2 48 , increment c = 1, and multiplier a = 3116728, 10430376854301, 
77596615844045 or 181465474592829. All these multipliers give good results on the spectral 
test in low dimensions, compared to other multipliers for the same modulus [43,44]. We 
verified that the runs with the four different multipliers gave results that are consistent 
within error bars for all the major observables. But when we analyzed the data for £/L, 
which ought to behave according to the finite-size scaling Ansatz 

£/L = x* + AL~ W + ... (4.2) 

where x* is a universal amplitude ratio and uj is a correction-to-scaling exponent, we found 
that our data fit this Ansatz very well (with uj = 0.82 from [39]) except for the points at 
L = 128 and L = 256, which showed deviations of magnitude 3% (~ 79 standard deviations) 
and 21% (« 170 standard deviations), respectively. Clearly something was going very wrong! 

After much work, we traced these systematic errors to the effects of long-range correla- 
tions (at lags that are multiples of large powers of 2) in the random- number generator [45-48]. 
It turns out [48] that these long-range correlations can arise within a single bond-update 
half-sweep of the Swendsen-Wang algorithm, provided that the lattice size is large enough 
compared to the modulus of the random-number generator. In a separate paper [48] we 
have studied these systematic errors in detail, in an effort to determine their approximate 
magnitude as a function of the lattice size and the random-number-generator modulus. Suf- 
fice it to say here that the systematic errors with a 48-bit random-number generator are 
comparable to or larger than our statistical errors only when the lattice size is a multiple 
of 64, which in this paper means L = 64, 128, 192, 256. We therefore discarded all the data 
for these lattices (~9.5 years CPU time, alas!) and performed new runs using 60-bit, 63-bit 
and 64-bit random-number generators: 
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Modulus m = 2 60 , multiplier a = 454339144066433781. 

Modulus m = 2 63 , multiplier a = 9219741426499971445. 

Modulus m = 2 64 , multiplier a = 3202034522624059733. 

(All these multipliers give good results on the spectral test in low dimensions, compared to 
other multipliers for the same modulus [43,44].) We also performed some additional runs on 
the smaller lattices using these generators. We have convinced ourselves [48] that generators 
using > 60 bits will exhibit significant systematic errors only on lattices larger than L = 256; 
in addition, they may exhibit slight systematic errors, less than about 2a, also at L = 256 
(we are currently investigating this latter issue more carefully). 

We assure the reader that the data reported in the present paper include only runs using 
"safe" random-number generators, i.e. m > 2 60 for L = 64, 128, 192, 256, and m > 2 48 for all 
other L. The CPU time figure of 17.8 years refers to these "good" runs only. 



5 Data analysis 

For each quantity O, we carry out a fit to the power-law Ansatz O = AL P using the 
standard weighted least-squares method. As a precaution against corrections to scaling, 
we impose a lower cutoff L > L m i n on the data points admitted in the fit, and we study 
systematically the effects of varying L min on the estimates of A and p and on the x 2 value. 
In general, our preferred fit corresponds to the smallest L min for which the goodness of fit is 
reasonable (e.g., the confidence level 14 is > 10-20%), and for which subsequent increases in 
£mm do not cause the \ 2 t° drop vastly more than one unit per degree of freedom. 

The behavior of the static quantities will be discussed in a separate paper [31]. Here we 
limit attention to the dynamic quantities. 

5.1 Integrated autocorrelation times 

Let us begin by summarizing the qualitative behavior of the integrated autocorrelation 
times T int) e> for different observables O, as reported in Tables El and EJ The three "energy-like" 
observables M, £, £' satisfy 

Tint,JV < nnt,E < Tint,£' (5-1) 

and the two "susceptibility-like" observables J\A 2 ,S 2 satisfy 

r mt,.M 2 < r int,5 2 j (5-2) 

in accordance with a rigorous theorem [11,25]. These five observables all have autocorrelation 
times in the same ballpark, as do C\ and Sq. The autocorrelation times of C2 and C3, by 

14 "Confidence level" is the probability that \ 2 would exceed the observed value, assuming that the under- 
lying statistical model is correct. An unusually low confidence level (e.g., less than 5%) thus suggests that 
the underlying statistical model is incorrect — the most likely cause of which would be corrections to scaling 
not included in the Ansatz. 
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contrast, are notably smaller. Of all the observables we measured, £' exhibits the largest 
autocorrelation time, with £ and S2 only slightly behind. 

Let us now fit the integrated autocorrelation times for all these observables to a simple 
power law r int q = AL Zint <° . We shall show the case of £' in detail; all the other observables 
behave similarly. 

In Figure ^ we have made a log-log plot of T- mt ^' versus L. (Please note that the error 
bars are significantly smaller than the plot symbols.) The plot shows notable curvature, i.e. 
there are fairly strong corrections to scaling, at least for L < 64. Consequently, the least- 
squares fits with L min < 64 all have enormous x 2 (confidence level < 0.05%), reflecting the 
fact that for L < 64 the corrections to scaling are many times our (very small) error bars. 
For L min > 96, by contrast, the x 2 values are good, reflecting the fact that in this regime 
the corrections to scaling are comparable to or smaller than our error bars. Our preferred 
fit corresponds to L min = 96, and yields 

Zint,£> = 0.4588 ±0.0047 (5.3) 

with x 2 — 0.352 (2 DF, level = 83.8%); here the error bar is one standard deviation (i.e. 
confidence level ~ 68%). 

A similar pattern is obtained for all the other observables, with the curvature always in 
the same direction. In all cases our preferred fit corresponds to L min = 96; the results of 
these fits are reported in Table El All the observables except C 2 and C 3 have exponents z int in 
the vicinity 0.45 ± 0.03. It is conceivable that the true values of these exponents are in fact 
exactly equal; we do not know whether the small differences between the estimates represent 
real differences or are merely the residual effects of corrections to scaling. 

It is worth noting that the rigorous inequality (J5.1|) implies 

Zint,Af < z int,£ < z int,£' , (5.4) 

while the estimates in Table |S] show the opposite behavior. This strongly suggests that in 
fact we have 

z int,N' = z int,£ = Zint,£' , (5-5) 

in accordance with the "almost-theorem" proven in [11, Section 2.2], and that the deviations 
in Table 03 result from corrections to scaling. Unfortunately, we don't know which of these 
estimates is closest to the true value; but we are inclined to trust more the estimate coming 
from the slowest of these modes, i.e. £'. We therefore give as our final estimate 

Ztotjsr = z int , £ = z int>£ , = 0.459 ± 0.005 ± 0.025 , (5.6) 

where the first error bar represents statistical error (68% confidence interval) and the second 
error bar represents possible systematic error due to the residual effects of corrections to 
scaling (68% subjective confidence interval). 

The susceptibility-like observables M 2 and S2 , by contrast, do show the correct inequal- 
ity arising from ()5.2j) . Comparing the two estimated exponents, and bearing in mind the 
"almost-theorem" that they should be equal, we give as our final estimate 

z- mt ,M* = ^int,5 2 = 0.443 ± 0.005 ± 0.030 . (5.7) 
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The estimates ()5.6|) and ()5.7|) are consistent with each other, as well as with the estimates 
for 2i nt) 5 and z ixit)Cl ■ This suggests that the true values of all these exponents might be exactly 
equal. Only the estimates for z intjC2 and z ixit)Cz are significantly lower than the others; and 
even here, it is conceivable that the discrepancy again arises from corrections to scaling. 



5.2 Exponential autocorrelation time and autocorrelation func- 
tions 

Recall that exponential autocorrelation time of an observable O is defined as 

-Itl , , 

T e xp,o = hm 7- , 5.8 

t^oo log poo (t) 

and that the exponential autocorrelation time of the system is defined as 

r exp = supr cxp . . (5.9) 
o 

All observables that are not orthogonal to the system's slowest mode satisfy r exPt o = r exp . 
Since all the observables studied in this paper are invariant under the symmetry group of 
the Potts model, we have no reason to expect that any of them are orthogonal to the slowest 
mode. We therefore expect — and will verify numerically — that they all have the same 
exponential autocorrelation time r exPi e>, which is presumably equal to r cxp . 

We shall begin by discussing the qualitative behavior of the autocorrelation functions for 
various observables at fixed L. Then we shall discuss the L-dependence of various quantities 
associated to the exponential decay of the autocorrelation functions. Finally, in the next 
subsection, we shall discuss the finite-size scaling of the autocorrelation functions. 

The typical behavior of the autocorrelation functions pooif) is depicted in Figure El For 
simplicity we have shown only the two observables exhibiting the most extreme behavior 
(among these we have measured): namely, £', which has the largest r int and whose auto- 
correlation function shows the least curvature (i.e. is closest to a pure exponential); and 
C2, which has the smallest T- mt and whose autocorrelation function shows the most curva- 
ture. The plots for all other observables are intermediate between these two. 15 Clearly, each 
autocorrelation function behaves asymptotically for large t as 

Poo{t) « A e- W/T ^° . (5.10) 

We obtained rough estimates of r cxPt o and the amplitude Ao by performing an unweighted 
least-squares fit to 

log Poo(t) = a-bt (5.11) 



15 These behaviors underlie our choice of the window factors for estimating r m t t o- namely, c = 8 for J\T, 
£, £', A4 2 , S2 and C\, whose autocorrelation functions are close to a pure exponential; c = 10 for So; and 
c = 15 for C2 and C3. 
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over the range 7"i n t,f < £ < 3Ti ntj £/ where all the autocorrelation functions are approximately 
a pure exponential; this yields r exp a = 1/6 and A = e aW The results are shown in Tables |H1 
and [7| Clearly, all the observables studied here have the same value of T exPi e>, as expected 
theoretically We shall use r cxp£ i from Table El as our best estimate of r exp . 

In Figure El we have plotted the estimated r exp versus L. We attempted to estimate the 
dynamic critical exponent z exp by fitting 

r exp w BL*~» . (5.12) 

In performing this fit, we used as rough error bars on r exp the sum of the error bar on Ti n t,£' 
and the standard deviation of the r cxp estimates for the seven observables M, £, £', Ai 2 , S 2 , 
S and d. Our preferred fit has L min = 96, and yields z exp = 0.481±0.007, B = 1.706±0.058 
(x 2 = 0.754, 2 DF, level = 68.6%). Of course, these error bars should not be taken terribly 
seriously. 

It is not clear whether z exp is equal to z int £> or is slightly larger. This question is related 
to the degree of curvature in the plot of the autocorrelation function, and more specifically 
to its L-dependence. To investigate this question in more detail, we first observed that 

^e-ltlA-p < Poo {t) < e - |t|/re * p (5.13) 

(as is obvious from Figure |2J). Therefore, if we define the modified autocorrelation time 

f m = - T e-l*l/— = - 1 + 6 u , (5.14) 
2 ^ 2 1 - e-V^p ' 

t=— oo 



we necessarily have 



Ao < ^ < i . (5.15) 

''"exp 

We therefore studied the L-dependence of the quantities Aq and Rq = T inti e> / ^cxp for various 
observables O (see Tables El and EJ) • 

We tried fits of Aq and Rq to the alternative Ansatze cL~ p and (with u = 0.82). 

Unfortunately, we do not have any valid error bars on Ao and Ro] but we can assign fictitious 
error bars and compare the relative \ 2 for the two fits. We did this for the two extreme 
observables, £' and Ci- In all cases we found that the power-law Ansatz gives a much better 
fit, and also one that holds over a wider range of L. We estimated p ~ 0.023 for Agi, 
p w 0.021 for R £ i, p m 0.092 for A c , 2 , p ~ 0.135 for R C2 . The values for p for A £ i and R £ > 
are very nearly equal, and are almost exactly equal to our estimate of z exp — Zmt,£' ~ 0.022. 
As for the values for p for Aq 2 and Rc 2 , they violate the rigorous inequality p(Ac 2 ) > p(Rc 2 ) 
that follows from ()5.15|) ; this strongly suggests that the two exponents p(Ac 2 ) and p(Rc 2 ) 
are in fact equal, though we do not know whether the correct value lies nearer to 0.092 or 
to 0.135. The latter value is fairly close to our estimate of z exp — ^int,e 2 ~ 0.127. 

16 In principle we should have done a properly weighted least-squares fit, taking account of the covariance 
matrix 1|3.10J) among the estimates poo(t)- But this is rather complicated, and we did not think it was worth 
the effort. As a result of this laziness, we are unable to quote any reliable error bars on T exP) or A. 
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Since for each of these observables it appears that Aq and Rq have the same exponent p, 
we tried fits of Rq/Aq to the Ansatz c\ + CiL~^ (with u = 0.82). For £' we find a limiting 
value Rs'/Agi ~ 1.010; for C 2 it is more difficult to tell, but a limiting value Rc 2 /A C2 ~ 1.14 
seems plausible. 



5.3 Finite-size scaling of autocorrelation functions 

A final way to study these questions is to investigate the finite-size scaling of the autocor- 
relation functions. The standard dynamic finite-size-scaling Ansatz for the autocorrelation 
function pooif) is 

P oo(t;L) « \t\^ho(—;^) ■ (5.16) 

(Here the dependence on the coupling constants, e.g. the inverse temperature, has been 
suppressed for notational simplicity.) Summing ()5.16j) over t, it follows that 

Tmt.o ~ r^ p P o , (5.17) 

or equivalently, 

Zint,0 = (1 -PoKxp.O ■ (5.18) 

Thus, only when p a = do we have z- int) o = ^ex P ,o [4] • In this latter case the Ansatz (|5.16j) 
can be rewritten in the equivalent form 



P oo(t;L) » h (—;^-) ■ (5.19) 

\7int,C> Li J 

To test this latter Ansatz, we have plotted log pooif) versus t/r inti o for the observable 
O — £' (Figure H]). For clarity we have included only the data from L > 12; the data coming 
from different lattice sizes are plotted with different symbols. We have also depicted for 
reference a line corresponding to the pure exponential p£>£> (t) = e~*/ Tint . £ ' . 

The data fall roughly onto a single curve, but there are clear corrections to scaling: the 
points move upwards (away from the pure exponential line) as L increases. It is not clear 
whether the points are tending to a limiting curve as L — ► oo, or whether they will continue 
indefinitely to move upwards. This is another way of saying that we do not know whether 
Pe> = 1 — Zint,£'/ Z exp is exactly zero or is slightly positive (e.g. ~ 0.067). 



6 Discussion 

In this paper we have obtained high-precision data at the critical point of the three- 
dimensional Ising model on fairly large lattices (up to L — 256), which have allowed us 
to derive quite accurate estimates of the dynamic critical exponents z intj0 and z exp for the 
Swendsen-Wang algorithm for this model. Our data resolve the discrepancies between pre- 
vious works [9,14-20], which can now be understood as arising from corrections to scaling. 

We would like to conclude by comparing our numerical results with some of the theo- 
retical frameworks that have been proposed by previous authors. These frameworks have 
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as their goal to understand the dynamic critical behavior of the Swendsen-Wang algorithm 
for the various ferromagnetic Potts models, and in particular to relate the dynamic critical 
exponent (s) ^sw to the static critical exponents for the same models. The three most impor- 
tant of these frameworks are the Li-Sokal proof [25] and its extensions [11,12]; the scaling 
Ansatz of Klein, Ray and Tamayo [21]; and the empirically based conjectures of Coddington 
and Baillie [16]. 

We have discussed the Li-Sokal bound zsw > a / v in the Introduction, and there is not 
much more to say. Suffice it to observe once again that while this bound is close to sharp 
(and possibly even sharp modulo a logarithm) for the two-dimensional Potts models with 
q = 2, 3, 4, it is clearly far from sharp in the three- and four- dimensional Ising models. Our 
data for the three-dimensional Ising model yield z$w ~ 0.46, compared to a/i/~ 0.1756 [41]; 
and for the four- dimensional Ising model it is generally believed that ^sw — 1 [16,21-24], 
compared to a/u = (x log 1 / 3 ). Clearly, some other physical mechanism, beyond the one 
exploited in the Li-Sokal proof, must be principally responsible for the critical slowing- 
down in these latter models; the central open problem is to identify this mechanism and to 
determine theoretically the dynamic critical exponent. 

Klein, Ray and Tamayo [21] have presented a scaling Ansatz leading to the conjecture 

zsw = z G - - — , (6.1) 

d m u 

where zq is the dynamic critical exponent for the Glauber dynamics in the same model, and 
d m is "the mean fractal dimension of the finite clusters" in the Fortuin-Kasteleyn represen- 
tation of the model. 17 The trouble with this Ansatz, alas, is not simply that the numerical 
value of d m is unknown; it is, rather, that the definition of d m is too vague to serve even 
as a guide for numerical attempts to determine its value. (The situation would be different 
if, for example, d m could be defined as a dimension associated with the scaling behavior of 
some specific observable.) Consequently, Klein, Ray and Tamayo were limited in practice 
to observing that d m presumably "lies between d, the spatial dimension, and df = d — fi/v, 
the fractal dimension of the incipient infinite cluster" [21, p. 164]. This plausible reasoning 
yields the conjectured inequality 

27 27 . , 

z g - ~j n < zsw < z G - — . (6.2) 

dv — p dv 

In Table [TO] we compare zsw with the Klein-Ray-Tamayo bounds %kx an d ^krt^ f° r 
the ferromagnetic Potts models (in dimension d < 4) having a second-order transition. The 
bounds appear to be violated for all three two-dimensional Potts models: for the Ising and 
3-state models, the violation is possibly within the errors in the determination of zsw an d 
z G , especially if one takes into account the possibility of logarithms; but for the 4-state 
model, the violation is blatant (barring a gross error in the determination of zq). For the 
three-dimensional Ising model, curiously, the lower bound seems to be exact (within errors); 

17 The discussion of Klein, Ray and Tamayo [21] does not distinguish between the various dynamic critical 
exponents Zi n t,0 and z cxp . For simplicity we can assume that their discussion refers to z 0X p for both the 
Swendsen-Wang and Glauber algorithms. 
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while for the four-dimensional Ising model, the upper bound is exact (modulo logarithms). 
It would be interesting to know whether the latter facts are anything more than curious 
coincidences. 

Coddington and Baillie [16] carried out a careful numerical study of the dynamic critical 
behavior of the Swendsen-Wang algorithm for the Ising models in dimensions d = 2,3, 4, on 
the basis of which they made the remarkable conjecture that for these models zsw = fl/v 
exactly. More specifically, they observed that the mean size of the largest cluster, here 
denoted C\, scales at the critical point as C\ ~ L d ~ l3 / U , and they found that their data could 
be explained by the asymptotic Ansatz 

T SW d/L d » a + MogL. (6.3) 

If true, this would imply that the correct dynamic critical exponent ^sw f° r the two- 
dimensional Ising model is neither (log) [64] nor 0.222 ±0.007 [10], but rather 1/8 (possibly 
multiplied by a logarithm). The data of [10] are certainly consistent with this possibility, 
though they do not distinguish it from the other possible behaviors. 18 

For the three-dimensional Ising model, the Coddington-Baillie conjecture would imply 
that zsw = 0.5183(4), which at first sight is incompatible at the 3a level with our estimate 
z mt,,£' — 0.459 ± 0.005 ± 0.025 (central value ± statistical error ± systematic error). Indeed, 
the curvature in Figure ^ assuming that it continues in the same direction, suggests that the 
true z- m t,£' is, if anything, slightly lower than our estimate based on L < 256. On the other 
hand, the difference between these exponents is small, and a small exponent is very difficult to 
distinguish from zero. We therefore made a direct test of the Coddington-Baillie conjecture 
by studying the combination T mt ^iC\/ L d . (Since we don't have statistically valid error bars 
for this combination, we used the triangle inequality to set worst-case error bars.) A fit to 
T mt,£'Ci/ L d = AL P yields a decent x 2 OT ^Y f° r L min > 96; our preferred fit is L min = 96 and 
yields p = -0.0573 ± 0.0052, log A = 0.670 ± 0.025 (x 2 = 0.261, 2 DF, level = 87.8%). This 
estimate for p is, not surprisingly, in almost perfect agreement with the values zsw = 0.459 
and p/u = 0.5183. 19 On the other hand, if we fit to T mt:£ ,Ci/L d = A + BL~ U with tu = 0.82, 
a decent x 2 is again obtained only for L min > 96; our preferred fit is again L min = 96 and we 
get A = 1.362 ±0.011, B = 6.064 ±0.553 (x 2 = 1-726, 2 DF, level = 42.2%). Figure shows 
the data points (□) and the corresponding fit. The fact that a reasonable fit is obtained 
with A far from zero (on the scale set by the observed values of T^giCi/L ) means that our 
data are also consistent with a behavior T- mt ^rCi/ L d — > A > as L — > oo, and hence p = 0. 

It is very hard (if not impossible) to distinguish, on purely numerical grounds, between 
these two behaviors. The x 2 is slightly better for the power-law fit, but this minor difference 
should not be taken terribly seriously. The bottom line, it seems to us, is this: the data 
shown in Figure El do not give any strong reason to believe that r- int £iCi/ L d is tending to 
zero as L — > oo. Indeed, inspection of the curve would suggest a limit in the range 1.2-1.3, 

18 Unfortunately, C\ was not measured in [10]. 

19 Just to be safe, we also checked the theoretically predicted scaling of C\. Using the Ansatz C\/L d = AL P , 
we get a good fit already with L min = 48: p = -0.5162 ± 0.0002, log A = 0.0892 ± 0.0008 (x 2 = 5.848, 4 DF, 
level = 21.1%). This value is not, strictly speaking, consistent with the estimate f3/v — 0.5183(4), but it is 
very close; the difference 0.002 can surely be understood as an effect of the residual corrections to scaling. 
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depending on the extent to which the curvature continues at larger L; this predicted limit 
is only about 20% below the maximum value attained at L 40, and is thus very far from 
zero. A more reliable judgment on the limiting value of r int ^rCi/ L d will have to wait 5-10 
years, when data will hopefully be available at (say) L = 512 and L = 1024. But as things 
stand today, our data are fully consistent with the Coddington-Baillie conjecture (albeit 
without a logarithmic term). 

If the Coddington-Baillie conjecture is interpreted as applying to z cxp rather than to 
Zint,£'> then the consistency between our data and the conjecture is even stronger. This is 
to be expected, as our estimate z exp ~ 0.481 is closer to the value (3/v ~ 0.5183. The data 
points for T exPj e>Ci/ L d are also shown (alas, without error bars) in Figure El (points *). The 
curvature is slightly weaker than for T^t£'C\/ 'L , and the data seem to be tending to a limit 
in the range 1.35-1.45. 

Of course, as Coddington and Baillie [16] themselves observe, z sw — fi/v cannot pos- 
sibly be a general identity for the Swendsen-Wang dynamics, as it clearly fails for the 2- 
dimensional Potts models with q = 3 and q = 4 (see Table ITU)) . Indeed, the Li-Sokal bound 
(II. 2|) ensures that we must have zsw > (3/v in any Potts model where a/v > (3/v. At best, 
the identity zsw = (3/v could hold for the special case of the Ising models. 

Since the Swendsen-Wang algorithm is defined naturally for all ferromagnetic Potts mod- 
els, a theoretical framework that is valid only for the Ising case seems unnatural and, in our 
opinion, unlikely to be correct. But the Coddington-Baillie conjecture can be rephrased in 
the following way so as to be potentially valid for all Potts ferromagnets. Suppose that there 
exists an as- yet-not-understood physical mechanism causing slowness of the Swendsen-Wang 
dynamics that is somehow related to the typical size of the largest cluster. In this case, an 
inequality of the form 

L d 

r S w > const x — , (6.4) 
Ci 

analogous to the Li-Sokal bound (jl.lj) . might hold for all Potts ferromagnets, irrespective 
of dimension and number of states. Indeed, it might even be possible to prove such an 
inequality rigorously (for one or another of the various autocorrelation times), if the physical 
basis were sufficiently well understood. Furthermore, it is even conceivable that the Li-Sokal 
mechanism and this new mechanism might together exhaust the reasons for slowness in the 
Swendsen-Wang dynamics, leading to the exact relation 

z S w — max(a/i/, (3/v) (6.5) 

(possibly modulo a logarithm) for all Potts ferromagnets. All currently available numerical 
data are consistent with the validity of the grand conjecture (J6.5|) . provided that a multi- 
plicative logarithm is permitted but not mandatory. 

This conjecture is, of course, a wild speculation; indeed, we consider it unlikely, a priori, 
for a dynamic critical exponent of any nontrivial dynamics to be exactly expressible in terms 
of static critical exponents (except for trivial cases such as Gaussian models). But stranger 
things have happened; and this conjecture is, in any case, certainly worth closer investigation. 
More modestly, this line of reasoning suggests that efforts be made to prove the inequality 
(16. 4|) and to understand what kind of physical mechanism might cause it to hold. 
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L 


Total # 


Total # 


Lengths (and numbers) 


T"int,f 


# measurements 


CPU time 


CPU time 




iterations 


discarded 


of individual runs 


in units of r int) £/ 


(/usec/spin) 


(years) 


4 


10 s 


10 5 


10 8 (1) 


2.3873 


4.2xl0 Y 


0.90 


0.000 


6 


10 s 


10 5 


10 s (l) 


3.1054 


3.2 x 10 7 


0.56 


0.000 


8 


10 8 


10 5 


lO 8 (l) 


3.7182 


2.7xl0 7 


0.45 


0.001 


12 


10 5 


10 5 


10 s (1) 


4.7469 


2.1 x 10 7 


0.40 


0.002 


16 


5xl0 8 


5xl0 5 


10 8 (5) 


5.6180 


8.9 xlO 7 


0.39 


0.025 


24 


10 s 


10 5 


10 8 (1) 


7.0712 


1.4xl0 7 


0.40 


0.018 


32 


2xl0 8 


2xl0 5 


10 8 (2) 


8.2563 


2.4 xlO 7 


0.42 


0.087 


A O 

48 


3.8 x 10 


2.1 x 10 


10 5 (2), 10' (17), 


10.2101 


3.7x 10 


0.43 


0.573 








5 xlO 6 (2) 










64 


3.1 xlO 8 


2.7xl0 6 


2xl0 7 (12), 10 7 (5), 


11.7976 


2.6 xlO 7 


0.44 


1.133 








2 xlO 6 (10) 










96 


8.6 xlO 7 


1.9 xlO 6 


5 xlO 6 (17), 5 xlO 5 (2) 


14.4677 


5.8 xlO 6 


0.45 


1.085 


128 


4.8 xlO 7 


2.4 xlO 6 


2.5 xlO 6 (16), 10 6 (8) 


16.5379 


2.8 xlO 6 


0.46 


1.467 


192 


2.9 xlO 7 


3.3 xlO 6 


10 6 (25), 5 xlO 5 (8) 


19.8511 


1.3 xlO 6 


0.47 


3.057 


256 


4.05 xlO 7 


4.6 xlO 6 


10 6 (35), 5 xlO 5 (11) 


22.6743 


1.6 xlO 6 


0.48 


10.335 



Table 1: Summary of our runs. The total CPU time used in these runs was approximately 
17.8 years. 



L 


X 


C H 


e 


E 


Ci/L a 


4 


21.1944 ± 


0.0033 


6.6297 ±0.0012 


2.52000 ±0.00036 


0.4310438 ±0.0000405 


0.512987 ±0.000043 


6 


49.0575 ± 


0.0089 


8.3079 ±0.0019 


3.81458 ±0.00057 


0.3881242 ±0.0000280 


0.424582 ± 0.000043 


8 


87.8535 ± 


0.0175 


9.4539 ±0.0024 


5.10153 ±0.00080 


0.3690531 ± 0.0000212 


0.368769 ± 0.000042 


12 


197.8047 ± 


0.0449 


11.0562 ±0.0032 


7.66909 ±0.00130 


0.3522118 ±0.0000140 


0.300979 ±0.000039 


16 


350.5792 ± 


0.0388 


12.2219 ±0.0017 


10.23802 ±0.00083 


0.3448934 ±0.0000047 


0.260170 ±0.000017 


24 


783.1500 ± 


0.2170 


13.9336 ±0.0050 


15.38221 ±0.00302 


0.3385100 ±0.0000068 


0.211590 ±0.000034 


32 


1382.5072 ± 


0.2921 


15.1928 ±0.0042 


20.52010 ±0.00303 


0.3357384 ±0.0000035 


0.182575 ±0.000023 


48 


3076.1461 ± 


0.5232 


17.0744 ±0.0038 


30.80276 ±0.00359 


0.3333254 ±0.0000016 


0.148217 ±0.000015 


64 


5420.9953 ± 


1.0953 


18.4697 ±0.0049 


41.08335 ±0.00563 


0.3322826 ±0.0000013 


0.127789 ±0.000015 


96 


12038.7855 ± 


5.1149 


20.5499 ±0.0114 


61.65067 ±0.01750 


0.3313744 ±0.0000016 


0.103645 ±0.000026 


128 


21193.9763 ± 13.0416 


22.1551 ±0.0179 


82.20142 ±0.03354 


0.3309822 ±0.0000016 


0.089301 ± 0.000033 


192 


47036.4986 ±41.9075 


24.4737 ±0.0288 


123.33211 ±0.07213 


0.3306421 ± 0.0000013 


0.072421 ±0.000039 


256 


83001.9797 ±66.8752 


26.2767 ±0.0280 


164.74988 ± 0.08662 


0.3304971 ± 0.0000008 


0.062502 ±0.000031 



Table 2: Static data from the Monte Carlo simulations at the critical point of the 3- 
dimensional Ising model. For each lattice size (L), we report the susceptibility (x), the 
specific heat (Ch), the second-moment correlation length (£), the energy (E), and the mean 
size of the largest cluster (C±). The quoted error bar corresponds to one standard deviation 
(i.e. confidence level w 68%). 
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L 




7int,£ 


T int,£' 


r int,A4 2 




A 

4 


oil «n -i- n nm q 
Z.llOy ± U.UUlo 


o Qf;o7 -i- n nnoi 
z.ooy < ± U.UUzl 


o qnnfi 4- n nnoo 
z.oyUo ± U.UUzZ 


o Q/inQ 4- n nnoi 
z.o4yo ± U.UUzl 


o QfiQn 4- n nnoo 
z.oooU ± U.UUzz 


P. 



o 7oc;7 -U n nnoft 

Z. ( 20 I ± U.UUZO 


q nfii s -U n nnQi 
O.UOlo ± U.UUOl 


q 1 nofi 4- n nnQi 
O.lUyo ± U.UUOl 


q nQm 4- n nnQi 
O.UOUl ± U.UUOl 


q nn77 4- n nnQ 1 
o.uy 1 1 ± U.UUOl 


o 
o 


q ooofi -U n nnQQ 
O.zzyo ± U.UUOO 


o.04yo ± U.UU4U 


q 70QS 4- n nn/i i 
0. ( zOo ± U.UU41 


q ^oc;o 4- n nnQn 
o.oyoy ± u.uuoy 


q 7nQn 4- n nn/i 1 
O. ( UOU ± U.UU41 


1 
LA 


4.U008 ± U.UU4< 


4.0014 ± U.UUOo 


a 7k^q 4- n nnfin 
4. < OOo ± U.UUOU 


a K'iOA 4- n nnc;fi 
4.00Z4 ± U.UUOO 


a 7i Q7 4- n nnc;n 
4. i is i ± u.uuoy 


1 P. 
10 


/i 77ni 4- n nno7 
4. 1 I Ul ± U.UUz / 


0.4ooo ± U.UUoo 


k fiQnQ 4- n nnQ/i 
O.OOUO ± U.UU04 


k Qnfin 4- n nnQi 
o.ouoy ± u.uuoi 


k c;7Qn 4- n nnQ/i 
o.o i oy ± u.uuo4 


0/1 
Z4 


oc;fi7 -U n nnsQ 
o.yoo i ± u.uuoo 


r fi/ins 4- n m no 
o.o4Uo ± U.UlUz 


7 nso7 xn m ns 
< .Uoy < ± U.UlUo 


p. c;7sn 4- n nnoK 
o.o ioi)± u.uuyo 


nnoQ 4- n m r\K 
o.yyzo ± u.uiuo 


32 


6.9303 ±0.0074 


7.9625 ±0.0090 


8.2727 ±0.0096 


7.5922 ±0.0084 


8. 1362 ±0.0094 


48 


8.5612 ±0.0073 


9.8308 ±0.0090 


10.2429 ±0.0096 


9.2535 ±0.0083 


10.0271 ±0.0093 


64 


9.8955 ±0.0101 


11.3374 ±0.0124 


11.8272 ±0.0132 


10.5715 ± 0.0112 


11.5342 ±0.0127 


96 


12.1713 ±0.0263 


13.8983 ±0.0321 


14.5123 ±0.0343 


12.7790 ±0.0284 


14.0740 ±0.0327 


128 


13.9684 ±0.0438 


15.8963 ±0.0533 


16.5963 ±0.0567 


14.4798 ±0.0462 


16.0373 ±0.0540 


192 


16.8861 ±0.0778 


19.0981 ±0.0933 


19.9312 ±0.0996 


17.0821 ±0.0790 


19.0814 ±0.0933 


256 


19.3990 ±0.0813 


21.8291 ±0.0969 


22.7669 ±0.1034 


19.4148 ±0.0814 


21.7656 ±0.0967 



Table 3: Autocorrelation times for the Swendsen-Wang algorithm at the critical point of the 
3-dimensional Ising model. For each lattice size (L) , we report the integrated autocorrelation 
time for the bond occupation (r int ^), the energy (r int ^), the nearest-neighbor connectivity 
( r mt,£'), the squared magnetization (7^x2), and the mean-square cluster size (Ti n t,s 2 ). The 
quoted error bar corresponds to one standard deviation (i.e. confidence level ~ 68%). 



L 


Tint ,<So 


TintA 


T mt,C 2 


r int,C 3 


4 


2.0347 ±0.0019 


2.2994 ±0.0020 


1.1843 ±0.0010 


1.4474 ±0.0014 


6 


2.5473 ±0.0026 


2.9784 ±0.0029 


1.3890 ±0.0013 


1.7324 ±0.0018 


8 


2.9917 ±0.0033 


3.5488 ±0.0039 


1.5639 ±0.0015 


1.9720 ±0.0022 


12 


3.7408 ±0.0046 


4.5025 ±0.0055 


1.8564 ±0.0020 


2.3786 ±0.0029 


16 


4.3847 ±0.0026 


5.3001 ±0.0031 


2.0973 ±0.0011 


2.7099 ±0.0016 


24 


5.4764 ±0.0082 


6.6205 ±0.0097 


2.4807 ±0.0031 


3.2554 ±0.0046 


32 


6.3760 ±0.0072 


7.6876 ±0.0086 


2.7875 ±0.0026 


3.6901 ±0.0039 


48 


7.8999 ±0.0072 


9.4380 ±0.0085 


3.2661 ±0.0024 


4.3843 ±0.0037 


64 


9.1455 ±0.0100 


10.8384±0.0116 


3.6451 ±0.0031 


4.9386 ±0.0049 


96 


11.2938 ±0.0262 


13.1816 ±0.0297 


4.2519 ±0.0074 


5.8452 ±0.0120 


128 


12.9878 ±0.0439 


14.9924 ±0.0487 


4.7240 ±0.0118 


6.5409 ±0.0193 


192 


15.7728 ±0.0783 


17.8138 ±0.0842 


5.4151 ±0.0194 


7.5908 ±0.0320 


256 


18.1480 ±0.0822 


20.2905 ±0.0870 


6.0248 ±0.0193 


8.5341 ±0.0324 



Table 4: Autocorrelation times for the Swendsen-Wang algorithm at the critical point of the 
3-dimensional Ising model. For each lattice size (L) , we report the integrated autocorrelation 
time for the number of clusters (T int < 5 ) and for the sizes of the three largest clusters (T int c- 
for i = 1,2,3). The quoted error bar corresponds to one standard deviation (i.e. confidence 
level « 68%). 
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Exponent 


Estimate 




X 2 (DF, CL) 


z 'mt,Af 

Zmt,S 

Zmt,£' 


0.4745 ± 0.0044 
0.4599 ± 0.0046 
0.4588 ± 0.0047 


96 
96 
96 


0.263 (2 DF, 87.7%) 
0.338 (2 DF, 84.4%) 
0.352 (2 DF, 83.8%) 


Zint,M 2 


0.4245 ± 0.0044 
0.4432 ± 0.0046 


96 
96 


1.641 (2 DF, 44.0%) 
1.126 (2 DF, 57.0%) 


Zmt,d 
z mt,C 2 
z int,C 3 


0.4832 ± 0.0047 
0.4384 ± 0.0045 
0.3537 ± 0.0034 
0.3836 ± 0.0040 


96 
96 
96 
96 


0.082 (2 DF, 96.0%) 
0.971 (2 DF, 61.5%) 
2.786 (2 DF, 24.8%) 
1.975 (2 DF, 37.2%) 



Table 5: Numerical estimates for the dynamic critical exponents of the 3-dimensional Ising 
model, based on least-squares fits with the specified L min . The quoted error bar corresponds 
to one standard deviation (i.e. confidence level 68%). 
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L 


Am 


T~exp,M 




Tcxp,£ 


A £ , 


7~exp,£' 


Am 2 


T exp,M 2 


As 2 


"?~exp,S2 


4 


0.8496 


2.3588 


0.9843 


2.3628 


0.9975 


2.3611 


0.9738 


2.3604 


0.9937 


2.3602 


6 


0.8372 


3.1071 


0.9704 


3.1059 


0.9895 


3.1080 


0.9561 


3.1063 


0.9845 


3.1078 


8 


0.8284 


3.7396 


0.9622 


3.7414 


0.9875 


3.7407 


0.9440 


3.7411 


0.9810 


3.7397 


12 


0.8157 


4.8034 


0.9524 


4.8069 


0.9844 


4.8060 


0.9275 


4.8062 


0.9748 


4.8065 


16 


0.8028 


5.7303 


0.9406 


5.7284 


0.9755 


5.7308 


0.9090 


5.7305 


0.9638 


5.7323 


24 


0.7919 


7.2722 


0.9283 


7.2692 


0.9682 


7.2687 


0.8899 


7.2575 


0.9548 


7.2621 


32 


0.7857 


8.5472 


0.9196 


8.5443 


0.9615 


8.5429 


0.8728 


8.5374 


0.9451 


8.5379 


48 


0.7827 


10.6314 


0.9112 


10.6380 


0.9548 


10.6364 


0.8532 


10.6385 


0.9333 


10.6400 


64 


0.7785 


12.3751 


0.9037 


12.3777 


0.9473 


12.3835 


0.8365 


12.4059 


0.9218 


12.4023 


96 


0.7709 


15.3193 


0.8879 


15.3373 


0.9326 


15.3276 


0.8126 


15.3481 


0.9041 


15.3353 


128 


0.7748 


17.6210 


0.8884 


17.6409 


0.9316 


17.6424 


0.8043 


17.6836 


0.8987 


17.6729 


192 


0.7782 


21.2902 


0.8867 


21.2979 


0.9289 


21.3083 


0.7949 


21.2413 


0.8933 


21.2720 


256 


0.7679 


24.6645 


0.8725 


24.6212 


0.9133 


24.6273 


0.7728 


24.6522 


0.8729 


24.6519 



Table 6: Estimates of the amplitude A and the exponential autocorrelation time T cxPi e> for 
the observables M , 8, £', M 2 and <S 2 - 



L 


A Sn 


"?~exp,.So 


Ac, 


7exp,Ci 


Ac 2 


r cxp, C'2 


Ac, 


T~exp,C:j 


4 


0.7991 


2.3582 


0.9479 


2.3499 


0.3298 


2.3284 


0.4858 


2.3399 


6 


0.7631 


3.1055 


0.9342 


3.0971 


0.3028 


3.0851 


0.4494 


3.0761 


8 


0.7500 


3.7358 


0.9297 


3.7252 


0.2931 


3.7079 


0.4318 


3.7036 


12 


0.7361 


4.7983 


0.9207 


4.7843 


0.2781 


4.7952 


0.4108 


4.7733 


16 


0.7248 


5.7269 


0.9078 


5.7042 


0.2698 


5.7093 


0.3970 


5.6844 


24 


0.7165 


7.2712 


0.8936 


7.2358 


0.2612 


7.2010 


0.3809 


7.2192 


32 


0.7131 


8.5471 


0.8838 


8.4986 


0.2520 


8.5026 


0.3710 


8.4798 


48 


0.7139 


10.6296 


0.8704 


10.5925 


0.2425 


10.6074 


0.3582 


10.5656 


64 


0.7125 


12.3728 


0.8585 


12.3428 


0.2353 


12.3740 


0.3471 


12.3565 


96 


0.7085 


15.3250 


0.8406 


15.2511 


0.2250 


15.3425 


0.3358 


15.2290 


128 


0.7154 


17.6157 


0.8335 


17.5967 


0.2214 


17.6720 


0.3301 


17.5485 


192 


0.7218 


21.2930 


0.8265 


21.2023 


0.2146 


21.3083 


0.3242 


21.0661 


256 


0.7139 


24.6895 


0.8078 


24.5465 


0.2076 


24.7020 


0.3105 


24.6323 



Table 7: Estimates of the amplitude Aq and the exponential autocorrelation time T ex p,e>. 
for the observables S , C\, C 2 and C 3 . 
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L 


Am 


Rm 


A £ 


Re 






Am 2 


R M 2 


A S2 


R S2 


4 


0.8496 


0.8834 


0.9843 


0.9889 


0.9975 


0.9976 


0.9738 


0.9804 


0.9937 


0.9945 


6 


0.8372 


0.8695 


0.9704 


0.9767 


0.9895 


0.9920 


0.9561 


0.9666 


0.9845 


0.9882 


8 


0.8284 


0.8583 


0.9622 


0.9699 


0.9875 


0.9896 


0.9440 


0.9556 


0.9810 


0.9841 


12 


0.8157 


0.8425 


0.9524 


0.9602 


0.9844 


0.9860 


0.9275 


0.9397 


0.9748 


0.9783 


16 


0.8028 


0.8303 


0.9406 


0.9501 


0.9755 


0.9800 


0.9090 


0.9237 


0.9638 


0.9702 


24 


0.7919 


0.8182 


0.9283 


0.9397 


0.9682 


0.9738 


0.8899 


0.9037 


0.9548 


0.9605 


32 


0.7857 


0.8103 


0.9196 


0.9310 


0.9615 


0.9673 


0.8728 


0.8877 


0.9451 


0.9513 


48 


0.7827 


0.8043 


0.9112 


0.9236 


0.9548 


0.9623 


0.8532 


0.8693 


0.9333 


0.9420 


64 


0.7785 


0.7987 


0.9037 


0.9150 


0.9473 


0.9546 


0.8365 


0.8532 


0.9218 


0.9309 


96 


0.7709 


0.7938 


0.8879 


0.9064 


0.9326 


0.9465 


0.8126 


0.8334 


0.9041 


0.9179 


128 


0.7748 


0.7915 


0.8884 


0.9008 


0.9316 


0.9405 


0.8043 


0.8205 


0.8987 


0.9088 


192 


0.7782 


0.7923 


0.8867 


0.8961 


0.9289 


0.9352 


0.7949 


0.8015 


0.8933 


0.8953 


256 


0.7679 


0.7876 


0.8725 


0.8863 


0.9133 


0.9243 


0.7728 


0.7882 


0.8729 


0.8837 



Table 8: Estimates of the amplitude A a and the ratio R = T int)0 /r cxp for the observables 
A/", £, £', M 2 and S 2 . 



L 


A So 


Rs 


Ac, 




A C2 


Rc 2 


Ac, 


Rc 3 


4 


0.7991 


0.8491 


0.9479 


0.9596 


0.3298 


0.4942 


0.4858 


0.6040 


6 


0.7631 


0.8126 


0.9342 


0.9501 


0.3028 


0.4431 


0.4494 


0.5527 


8 


0.7500 


0.7950 


0.9297 


0.9431 


0.2931 


0.4156 


0.4318 


0.5241 


12 


0.7361 


0.7756 


0.9207 


0.9335 


0.2781 


0.3849 


0.4108 


0.4931 


16 


0.7248 


0.7632 


0.9078 


0.9225 


0.2697 


0.3650 


0.3970 


0.4717 


24 


0.7165 


0.7522 


0.8936 


0.9094 


0.2612 


0.3407 


0.3809 


0.4472 


32 


0.7131 


0.7455 


0.8838 


0.8989 


0.2520 


0.3259 


0.3710 


0.4315 


48 


0.7139 


0.7422 


0.8704 


0.8867 


0.2425 


0.3068 


0.3582 


0.4119 


64 


0.7125 


0.7381 


0.8585 


0.8748 


0.2353 


0.2942 


0.3471 


0.3986 


96 


0.7085 


0.7366 


0.8406 


0.8597 


0.2250 


0.2773 


0.3358 


0.3812 


128 


0.7154 


0.7360 


0.8335 


0.8496 


0.2214 


0.2677 


0.3301 


0.3706 


192 


0.7218 


0.7401 


0.8265 


0.8358 


0.2146 


0.2541 


0.3242 


0.3562 


256 


0.7139 


0.7368 


0.8078 


0.8238 


0.2076 


0.2446 


0.3105 


0.3465 



Table 9: Estimates of the amplitude A and the ratio Rq = T int)0 /T cxp for the observables 
So, Ci, C 2 and C 3 . 
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Model 




P/v 






(lower) 
%RT 




(upper) 
KRT 


d = 1, all q 







1 


2 











d = 2 Ising 


(xlog) 


1/8 


7/4 


2.16(2) 


0.29(2) 


0.222(7) [10] 


0.41(2) 


d = 2, g = 3 


2/5 


2/15 


26/15 


2.19(2) 


0.33(2) 


0.514(6) [11] 


0.46(2) 


d = 2,q = 4 


1 (xlog" 3/2 ) 


1/8 (xlog 1 / 16 ) 


7/4 (xlog" 1 / 8 ) 


2.25(6) 


0.38(6) 


1 (xlog'") [13] 


0.50(6) 


d = 3 Ising 


0.1756(25) 


0.5183(4) 


1.9634(8) 


2.04(4) 


0.46(4) 


0.46(3) [this work] 


0.73(4) 


d = 4 Ising 


(xlog 1 / 3 ) 


1 (xlog 1 / 6 ) 


2 


2 


2/3 


1 (x log YY ) 


1 



Table 10: Dynamic critical exponent z- mt ,s' of the Swendsen-Wang algorithm for various 
Potts models, compared to the exponents arising in the theoretical frameworks of Li and 
Sokal [25], Klein, Ray and Tamayo [21], and Coddington and Baillie [16]. Static critical 
exponents a/is, (3/v and 7/z/ are exact values for the d — 1 and d = 2 models [13,28-30] 
and for d = 4 Ising [49], and are the best currently available numerical estimates for d = 3 
Ising [41]. Exponent zq for Glauber dynamics is taken from [50, Table 1] for the d = 2 
models (see also [51-55]), from [51,53,56-62] for d = 3 Ising, and from [63] for d = 4 Ising. 
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Least — squares fit to r 
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Autocorrelation functions for L=64 




20 40 60 80 fOO 

t 

Figure 2: Autocorrelation functions poo(t) for observables O — £' (□) and C 2 (*) at L — 64. 
We have also depicted the straight-line fits corresponding to r cxp C >. 



36 



Least — squares fit to r 
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Finite — size scaling of autocorrelation function 




2 4 6 8 10 

^/ T int,E' 



Figure 4: pe>£>(t) versus t/r hlt ^i for 12 < L < 256. The different symbols denote the different 
lattice sizes: L = 12 (+), L — 16 (x), L = 24 (□), L = 32 (0), L = 48 (o), L = 64 (*), 
L = 96 (X), L = 128 (X), L = 192 (-<>■), L = 256 (*). For comparison, we have also depicted 
the line corresponding to the pure exponential P£>£>(t) = exp(— t/r^e')- 
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Figure 5: The Coddington-Baillie products T int ^rCi/ L d (□ with triangle-inequality error 
bars) and r cxp ^/Ci/ L d (* without error bars) versus L~ - 82 . The least-squares fit line is 
T\nt,S'Ci/L d = 1.362 + 6.064L- - 82 and is obtained for L min = 96. 
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